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Abstract 

This is the second in a series of papers that investigates the topological sus- 
ceptibility in the interacting instanton liquid model (IILM); it deals with the 
technical issues relating to the Monte Carlo simulations that are specific to 
finite temperature. The IILM reduces field theory to a molecular dynamics de- 
scription, and for 'physical' quark masses the system behaves like a strongly 
associating fluid. We will argue that this is a generic feature for very light Dirac 
quarks in a non-trivial background, described in the semi-classical approach. To 
get rid of unnecessary complications, we will present the ideas of biased Monte 
Carlo, and implement the transition probabilities, for a toy model. 



1. Introduction 

In a recent paper [2Cj , we have performed simulations in the interacting in- 
stanton liquid model (IILM) at zero temperature. In this and the next paper 



2l| we will address finite temperature simulations. Using the light 'physical' 
quark masses obtained in [2(3], we have found that the Monte Carlo techniques 
that worked well at zero temperature become rather inefficient at finite tem- 
perature. The reason lies in the formation of instanton-anti-instanton pairs. 
They provide the mechanism for chiral symmetry restoration within the IILM 



17| : because the quark wavefunctions become localised on those pairs, 
the Dirac eigenvalue spectrum p(A) has no support at zero virtuality; according 
to the Banks-Casher relation [2j, 

(qq) = -np(\ = 0) , (1) 

chiral symmetry is restored. This is in contrast to zero temperature where 
the quark wavefunctions are delocalised and the eigenvalue spectrum extents 
down to zero virtuality. The localisation of the Dirac zero modes at finite 
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temperature is potentially troublesome for ordinary Monte Carlo techniques 
because it induces short-ranged interactions between instantons. 

In another numerical study [l6| of the IILM at finite temperature, no tech- 
nical problems were encountered because the quark mass parameters were large 
enough for ordinary Monte Carlo to work well. But as the masses decrease, the 
interactions become stronger and random sampling starts to run into trouble. 
Furthermore, in the deconfined phase screening sets in: the temperature fluc- 
tuations obstruct the formation of coherent field configurations that exceed the 
screening length, so that instanton sizes are cutoff at p < 1/T. Since the in- 
teractions in the IILM follow from the overlap of the instanton profiles, smaller 
sizes lead to shorter-ranged forces between pairs. 

The strong attraction is a generic feature of the IILM at finite temperature. 
Normalising the low frequency determinant to the dilute gas, the quark effective 
interaction is of the form 



Depending on the topological charge, the determinant is over the subspace 
spanned by the Nj or N A zero modes: T 2 = T^T if Nj > N A , and T 2 = TT^ 
otherwise. The matrix M 2 is diagonal and of the same dimension as T 2 ; its 
non- vanishing elements are given by the squared quark mass. Crucially, it is an 
attractive interaction because the determinant is bounded from below by unity. 
It is clear from the above expression that the quark interaction becomes ever 
stronger as quark masses decrease. 

The feature of strong and short-ranged interactions is thus not confined to 
the specific case of the trivial holonomy calorons [7fl that we have been study- 
ing, but is a general characteristic of non-trivial backgrounds; in particular, it 
will play a role for the newly found non-trivial holonomy calorons 0, El and 
(l2| . The latter might be the correct degrees of freedom to explain the con- 
finement/deconfinement phase transition , the lack of which is a major 
shortcoming of the current IILM. 

Systems with strong and localised interactions have been investigated for a 
long time in chemical engineering and computational chemistry, and are known 
to present computational challenges. They run under the name of strongly 
associating, or ionic, fluids. The technical problems that Monte Carlo methods 
face are two-fold: 

• The small relative volume of attraction makes ordinary Monte Carlo up- 
dates miss them most of the time. 

• Once a pair is formed, ordinary Monte Carlo moves can get stuck in these 
configurations because of the large energy difference. 

From an algorithmic point of view, it means that the acceptance rates to reach 
or leave the regions of phase space corresponding to instanton-anti-instanton 
molecules are very low. This would not be a problem if we had an infinite 
amount of computer time, but, in practice, it leads to very long autocorrelation 
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times. This, in turn, leads to strong dependence on initial conditions and, 
in the extreme case of very poor mixing between correlated and uncorrelated 
states, the system might stay fixed in one of those 'phases' during the available 
computer time. For all practical purposes, we lose ergodicity because we end 
up with samples that are either devoid or dominated by pairs, none of which is 
a representative sample. 

The solution is to develop algorithms which explicitly sample the attraction 
centres and are able to break up pairs with high acceptance rates. To sample 
these volumes, we need to, first, construct these regions and, then, to define a 
measure over them. This is a non-trivial geometrical problem; it also depends 
very strongly on the specific problem. Considering the dipole character of the 
instanton interactions, this seems like a daunting task. 

In recent years, general purpose algorithms have been developed which do 
not rely on an accurate construction of the union of all attraction volumes but 
on the different ways that a given, simple interaction box can be reached, e.g. 
the Unbonding-Bonding algorithm [22| . In its original form, it has been given 
for a canonical ensemble. We will adapt it to our needs for grand canonical 
simulations. 

As argued above, the quark interactions will generically lead to strong in- 
teractions for small quark masses, and any non-trivial background will lead 
to screening effects at finite temperature. To avoid unnecessary complications 
and background dependent features, e.g. the colour-orientation dependence for 
calorons, we will use a toy model that mimics the interactions of (J2j> . 

In section [2] we will quickly review the idea behind Markov chain Monte 
Carlo. We will then discuss biased Monte Carlo in section 02 and set up the 
framework for dealing with strong and short-ranged interactions in the grand 
canonical ensemble. Finally, we will present a toy model in section @J and 
benchmark the biased moves from section [3] within that setting. 

2. Markov Chain Monte Carlo 

The phase space for a general statistical mechanical system will be very high- 
dimensional for a large number of particles. Usual integration rules over a grid 
are not well adapted to evaluate the defining partition function because, for such 
high dimensional integrals, the number of grid points would be overwhelmingly 
large. 

Markov chain Monte Carlo techniques are the method of choic^l because 
they can 'find' the relevant regions that dominate the partition function. It is 
achieved by preferentially sweeping those parts of parameter space that have a 
large measure. It can be shown [13j that p™- — > p° q , where Pij is the transition 
probability and p cq is the unique invariant equilibrium distribution 
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So-called ergodic theorems relate ensemble averages to time averages of paths 
through phase space, i.e. 



1 



JV 



(4) 




with a specific path given by the Markov chain {X„}. 

In our case, we know p eq , essentially the integrand of the partition function, 
and we want to construct pij that converge to this equilibrium distribution. 
Using the convergence theorems, this is a well posed problem if we demand that 
p cq is invariant for p« . It is not hard to see that Q is fulfilled if we impose the 
stronger condition of detailed balance, 



This deceptively simple looking equation lies at the heart of all Monte Carlo 
simulations. Note that the Pij's are not unique and that there is a large amount 
of freedom in choosing them. This redundancy can be used to accelerate con- 
vergence. 

If there are different paths that connect states i and j, p^j = J^a Pij' © 
can be met by imposing the stronger condition 



This is sometimes called super-detailed-balance [5{. If there are very many 
different paths, N 3> 1, the implementation can become rather tedious, and 
super-detailed-balance might be the only viable option. Even if there are a 
manageable number of paths, the book-keeping needed to relate them might be 
overwhelming. 

It is important to note that the transition probability is really the prod- 
uct of the proposal probability Vij and the acceptance probability Aij. This 
observation lies at the heart of biased Monte Carlo techniques: we can tweak 
the proposal probabilities to increase acceptance rates. Typically, this leads to 
asymmetric transition probabilities. The latter can, however, also be found in 
ordinary Monte Carlo: consider for instance a simple insertion/deletion step 
for grand canonical simulations; the proposal probability for an insertion cor- 
responds to the probability to place the particle within the simulation box, 
whereas for a deletion it gives the probability to choose a particle already in the 
box. Clearly, these will be different in general. 

We will focus on the Metropolis algorithm which defines the acceptance 
probability by 



eq eq 
Pi Pij = PjPji ■ 



(5) 



cq a cq a 

Pi Pij =PjPji- 



(6) 





It is not hard to see that it satisfies ((SJ). 
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3. Biased Monte Carlo 



But for the simplest statistical systems, ordinary Monte Carlo techniques can 
become rather inefficient, if not downright inadequate: inefficiency manifests 
itself by long autocorrelation times, and inadequacy stems from an inability to 
sample the relevant regions accurately. For finite computer time, and hence 
finite-sized samples, we basically lose ergodicity. 

By 'ordinary' Monte Carlo moves, we mean random sampling. A typical 
example is updating the position of a particle: it is straightforward to implement 
such a move by adding a random increment to the particle's current position. 
This random increment is typically drawn from a fixed volume with uniform 
measure. Thus, V%j is constant for all updates and cancels in ([7]). 

For strong and short-ranged interaction, such a random increment will be 
very inefficient. Either we choose it so small that we can sample the interaction 
region, in which case very many sweeps are needed to move through phase space; 
this leads to very long autocorrelation times. Or we choose a large enough 
increment, so as to sweep quickly through phase space, but thereby missing 
the interaction regions most of the time; the sample will most likely not be 
representative. Furthermore, the acceptance probability for random sampling is 
typically given by the ratio of two Boltzmann factors, and does only depend on 
the energy difference. Due to the strong interaction, the acceptance probability 
will be very low for moves that attempt to leave the attraction centres. The 
Markov chain can become trapped in these energy-dominated configurations, 
which leads again to a non-representative sample. 

The way out is to use importance sampling. It is designed to sample those 
parts of phase space that dominate the partition function. It might help to get 
a rough criterion for importance sampling. Remember that in our case Monte 
Carlo methods try to evaluate the integral given by the partition function. If the 
integrand is peaked so strongly that the exact integral is well approximated by 
the integral over some localised patches, the algorithm will need to preferentially 
sample these parts of phase space. Thus, random sampling will be inefficient if 

V < AVexp(-H(AV)) , (8) 

where V is the volume of the simulation box, and AV <C V is the small region 
where the interaction is very strong. 

As mentioned in the introduction, in chemical engineering and computa- 
tional chemistry the issues of low acceptance rates have long been known: they 
are very low in polymer physics simulations, due to conformational obstructions; 
or in simulations of ionic fluids, due to strong short-ranged interactions. The 
latter is of immediate interest to us because the IILM at finite temperature dis- 
plays the same characteristics. To treat such systems correctly, new algorithms 
based on biased Monte Carlo have been developed [l9[ . We can take advantage 
of these well tested techniques in IILM simulations at finite temperature. 

Recently, efficient and general purpose algorithms have been developed 0] 



22[. In these algorithms, the focus is not on an accurate construction of the 



union of all the interaction regions, which is a difficult, and problem-dependent, 
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geometric task, but on the individual interaction regions and all the possible 
routes that lead to the same final state. In [U, the algorithm is further simplified 
by using super-detailed-balance whereas the algorithm in [22J does not rely 
on this stronger condition, and was shown to converge faster. We use the latter 
scheme, the Unbonding-Bonding algorithm (UB). 

The UB algorithm starts by defining a bonding region. It does not necessarily 
need to be the exact physical bonding region, but a strong departure from it 
will not be very useful. A list is made of those instantons that are in at least one 
bonding region of an anti-instanton and the same is done for the anti-instantons. 
There are Nf (iVjf ) bonded instantons (anti-instantons). Nf(ii) is the number 
of anti-instantons that instanton ij is bonded to, and analogously for N§(ia)- 
In what follows, ii means instanton ij, but can also stand for the state the 
instanton ij is in; unprimed quantities are evaluated before the move, whereas 
primed ones denote the same quantity after the move. 

We will now focus on the UB for instantons; the case for anti-instantons is 
then obvious. The bonding move consists of choosing uniformly an instanton 
and an anti-instanton, and placing the instanton in the bonding region of the 
anti-instanton with flat measure. The unbonding move consists of choosing one 
of the bonded instantons and to place it randomly in the simulation box; again 
both steps are performed with a flat measure for the bonded instantons and the 
simulation box respectively. This leads to the following transition probabilities 

n^M) NtN a Vi A ' v ' 

V -> = WV> (10) 

where (i/, %a) is a bonded instanton-anti-instantion pair, and Vi A is the bonding 
region of anti-instanton ia- The UB algorithm now adds up all possible routes 
that lead to the same final state i', bonding and unbonding; super-detailed- 
balance is, however, used with respect to the unbiased displacement move. The 
forward and backward proposal probabilities are then given by 

E p?«x,i*) + ^&z> (i 2 ) 

iA 

with Sf = 1 if i is bonded and 8f — otherwise. We assume that bonding and 
unbonding moves have the same a-priori-probability, 1/2, which we omitted 
because it cancels out anyway. 

Since we want to perform grand canonical simulations, we also need insertion 
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and deletion moves. The unbiased moves are given by 

1 } n' i ,n i +i = y 7 (13) 

p %+w = nt+t- (14) 

The corresponding biased insertions and deletions will be constructed along the 
lines of the UB algorithm, either by placing the instanton ij into the bonding 
region of an anti-instanton, or by removing the bonded instanton ij. More pre- 
cisely, insertions consist of choosing an anti-instanton i^, placing the instanton 
ij uniformly within the bonding box Vi A and finally summing over all possible 
anti-instantons that could have been chosen to reach this same final state; to 
delete an instanton, we choose uniformly from the list of bonded instantons Nf . 
In formulas, 



V N U N I+ X- E J^J' (15) 



B 



In this case, it does not produce much overhead to combine the biased and 
unbiased moves. More importantly, we found that acceptance rates could be 
boosted upon combining biased and unbiased insertions/deletions. To add them 
up, we need to specify the relative weights, the a-priori-probability for biased 
updates p&. The full insertion/deletion proposal probabilities are 

V Nl , Nl+1 = Pb V h NliNl+1 + (1 - Pb )Vf I>N[+1 , (17) 
V Nx+ i, Nl = Pb V h Nl+hNl + (1 - Pb)Vf I+1<Nj . (18) 

To reiterate, a similar factor for the a-priori-probabilities of canonical moves was 
tacitly omitted before because we chose to follow the original implementation 
of the UB algorithm, which does not mix biased and unbiased particle updates: 
its use of super-detailed-balance implies that the a-priori-probabilities give an 
overall multiplicative factor that drops out. 

Other than good mixing between bonded and unbonded structures, cluster 
moves have been argued to be important to accelerate convergence towards the 



equilibrium distribution IJ]. For the IILM, we assume that instanton-anti- 
instanton pairs are dominating in this respect. We therefore augment our list of 
single particle moves by pair-displacements and pair-insertions and -deletions. 

Whereas we made sure to have a uniform distribution among the Nf bonded 
instantons for the UB-type moves discussed so far, regardless of whether they 
were bonded many times, we naturally demand uniformity in the number of 
pairs Np for the pair-moves. 

The pair-displacements consist of collective translations of the pair, by an 
increment 5 <E vc, and of internal displacements, S £ Vj, of one of the pair's con- 
stituents. The latter are rejected if the displaced instanton leaves the bonding 
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box. The proposal probabilities for these moves are given by 

V c = 4r— , (19) 

N P v c v ' 

^ 7 = T^-, (20) 
N P vi 

We use super-detailed-balance for the internal displacement because the accep- 
tance rates could not be boosted by including unbiased and/or UB moves; the 
additional overhead, required to go beyond super-detailed-balance, is thus not 
justified. 

The pair-insertions and -deletions, set up without recourse to super-detailed- 
balance, are built from biased and unbiased moves, whose probabilities are 
summed up in the end. A biased pair-insertion consists of placing, uniformly, 
either an instanton %i or an anti-instanton ia in the simulation box; the partner 
is then positioned randomly within the bonding box. Note that the probabili- 
ties for both possibilities are added up. The unbiased pair-insertion consists of 
placing, uniformly, both an instanton ii and an anti-instanton ia in the simu- 
lation box. For a biased pair-deletion we select an instanton and anti-instanton 
uniformly from the pairs Np, whereas for an unbiased pair-deletion we select 
an instanton and anti-instanton randomly from Ni and Na respectively. The 
final proposal probabilities are then given by 

^ , rP /ll 1 1 1 1 \ , .,11 

W-*% + (i-ittj^- (22) 

The factor Sf iiA ensures that the biased contribution to the proposal probability 
is added only if %i and Ia are paired, i.e. Sf^ = 1 if instanton ii and anti- 
instanton %a are paired and 8f iiA — otherwise. The a-priori-probabilities p' b 
can be chosen to tune the acceptance rates even further. 

Given these different proposal probabilities, together with the trivial case 
of unbiased displacements, we can compute the acceptance probability through 
©. 

Finally, we need to specify the increment volumes Vi and the bonding boxes 
Vi, and we have to decide upon the different a-priori-probabilities pb, p' b and pi. 
The latter give the probability to perform a specific update and have, so far, 
been omitted because they will always cancel in the acceptance probabilities. 
These choices will depend on the problem at hand, and some fine-tuning runs 
cannot be avoided to set these parameters. 

In practice, we have finite-sized samples, and the outcome of the simulations 
will depend on these parameters. The differences will vanish with increasing 
sample size; this observation provides a straightforward means to test whether 
the biased updates have been implemented correctly. For computer-intensive 
simulations, large samples are prohibitive, and it is important that the depen- 
dence is rather weak if the algorithms are to be useful. 
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Figure 1: The solid lines represent the interaction for an instanton-anti-instanton pair. The 
repulsion for like-charged pairs is represented by dashed lines. It is chosen so strong that 
clusters larger than simple pairs are strongly disfavoured. We see that as the temperature t 
is raised the attractive well deepens. In this graph a = 2. 

4. Toy Model 

Now that we have set up the update moves using biased Monte Carlo tech- 
niques, we will test their performance compared to random sampling. We will 
use a two-dimensional toy model that nevertheless will mimic the IILM. The 
advantage of the toy model is that it will be computationally very cheap. Also, 
we can adjust the parameters freely to get a rough estimate for those regions 
in parameter space where pair formation will be important and, hence, biased 
Monte Carlo essential. 

The interactions in the IILM have the following form 

s » = b ( 1 + i)-T^- w ' h ( 1 + e - Mj!2 S) 

The logarithmic repulsion in the pair separation R, given here in units of the 
inverse temperature T, is typical for the IILlvd. It combines with the rational 
function to produce the large separation decay. The instanton-anti-instanton 
pairs feel an additional attraction due to the quark wavefunction overlaps. Their 
R 2 dependence follows from 'strong' overlaps, whereas the exponential describes 



2 In the ratio ansatz, see 
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the large separation behaviour. For our purposes, the product of both will be 
a good enough approximation. The quark mass is given by m, the number of 
flavours by Nr. As for the IILM, Su and Si a have the same fall-off behaviour. 

In two dimensions, this isotropic interaction results in a rather dense en- 
semble because the repulsion is mild and the attractive wells are comparatively 
deep. More importantly, these interactions will favour large clusters of instan- 
tons and anti-instantons. Such a behaviour does not correspond to the situation 
that we expect from the IILM, namely the formation of rather isolated pairs. 
In the IILM, the interaction is orientation dependent, and large conglomerates 
of instantons cannot form because most of the relative orientations within the 
cluster would be repulsive. To achieve this same effect of pair formation with 
our toy model, we add an extra 'entropic' repulsion term to S?J t . 

The exponential decay for the quark overlaps starts for e~ 27rR t 2 wl;we have 
set t = T/m, i.e. the temperature in units of the quark mass. For separations 
smaller than Rq = \n.t/w, the attraction wells are rather deep, and we choose the 
repulsion to become strong in order to obstruct the formation of large clusters. 
The functional form of this 'entropic' repulsion is rather unimportant, and for 
simplicity we choose it to be given by 

Sff=(^f)\ (25) 

It turns out that choosing a = 2 — 3 is sufficient to prevent cluster formation. 
For dilute ensembles, the precise value is rather irrelevant. 

Now that we have fixed the interactions, see also Fig.[TJ we define our system 
by the following partition function 



jNt ]N A 

dxNl+NA h.h^. expi ~ Sint{x)h (26) 



Ni N A Ni,N A 

q \ s runt i ocnt i \ * rant , ocnt , \ oint /97\ 

i<3 i<3 il^A 

where the free parameters of the model are d, t and the simulation box V. In 
the full model, instantons have a size p, and quantum effects are such that small 
sizes are disfavoured, d(p) — p 13 with j3 > [18[ ; the free parameter d plays that 
role. 

Neglecting interactions, d determines the density of the ensemble, 

^=2d. (28) 

This is the dilute gas result. Once interactions are included, the system will be 
more dilute if the total interaction is repulsive and less dilute otherwise; this is 
adjusted by the parameter t which captures the combined effects of temperature 
and quark masses, see (f23|) . 

In order to use the biased Monte Carlo framework set up in section [3l we 
need to decide on the bonding box to use. Given that the interaction is isotropic, 
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Figure 2: For the biased Monte Carlo moves to work well, the bonding box should be adapted 
to the interactions. Isotropy suggests the best form to be in the shape of an annulus. The 
shaded region corresponds to the bonding box. A priori R mln and i?, ma x are free parameters 
that can be fine-tuned to achieve good mixing, i.e. low autocorrelation times. 



we will choose an annulus as bonding box, see Fig. [2] It is natural to fix the 
free parameters i? m in and i? m ax so as to achieve low autocorrelation times. 

We perform biased simulations for three sets of bonding boxes, and compare 
them with unbiased runs. The natural box size extends to i? m ax = Ro, the point 
at which the quark interaction starts to vanish exponentially. The corresponding 
lower edge i? m i n is determined by Sj^R^n) ~ ^(^ma). Given the relatively 
large values for f, we can solve this approximately by 

Rnun « exp ^-Is} n j(i? max ) - ~ IntJ . (29) 

We also choose a larger interval with i? max = 2R and a smaller one that samples 
predominantly the very strong interaction region. In the latter case, we choose 
R min and R max according to SfXi^max) ~ Sf}(R B1 i n ) « ±mm R Sy i J[- A typical 
setup is displayed in Fig. [3] 

Before we start the simulations, it is worth mentioning one more technical 
point, related to the choice for the equilibrium probability density. We claim 
that it is given by 

jNj jNa -s iILt 
p0q = ^_e (3Q) 

We want to put attention to the fact that in we neglect the Boltzmann- 
counting factors for indistinguishable particles. We would like to argue that they 
are an artifact. These terms are introduced so as to allow for an unconstrained 
integration region, xf S [0, L a ], instead of the complicated multi-dimensional 
region given by {x^\x1 < x% ■ ■ ■ < x a N .}, with a labelling a component. The 
result of the integrations is corrected by dividing by the factorials because 
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Figure 3: We choose three different sets for (i2Lf_, Umax)' A small interval that samples the 
strong interaction region, a 'natural' one that sample most of the attractive well and a larger 
box that extends beyond what would naturally look like the bonding region. 

U CT {a:!? < x a{2) ' " < x a(N )} ~ IP 5 L a ] Nj , where a is a permutation. There- 
fore, the Boltzmann-counting factors are not really part of the weight A 
different, but equivalent, point of view is to keep the factorials in the weight and 
to sum the equilibrium probability, with the factorials included, over all permu- 
tations iVj! that correspond to mere relabellings of the particles in the initial 
state [l^|. The permutations cancel off the factorials, and we are left again with 
(|5D|) . The transition probabilities are similarly summed over all initial and 
final state permutations, and the resulting factor of Ni\Nj\ trivially cancels out 
of . The simplest way to understand the need to sum over permutations, lack 
of which will violate detailed balance, is to realise that, taking indistinguisha- 
bility at face value, the states are really equivalence classes^. Since we actually 
use a specific representative but the result should be independent of it, we need 
to sum over all members of the given equivalence class. The upshot is that, for 
all practical purposes, we treat the different particles as distinguishable! 

The parameter t is set to the values given in Fig. [TJ The dilution parameter 
d is chosen small enough so that we end up with a rather dilute system because 
the biased Monte Carlo moves, set up in section [31 will not work well for dense 
ensembles. We are interested in how the system behaves if d is further decreased. 
We can envisage three different equilibrium states: either the ensemble will 
be dominated by pairs, or equilibrates in a mixture of paired and unpaired 



3 The group by which we factor out is the permutation group 
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Figure 4: Choosing a rather small parameter for t, i.e. a shallow attraction well, the biased 
(bottom) and unbiased (top) simulations have thermalised over the same time scale. The over- 
head that importance sampling introduces is not necessary provided that the autocorrelation 
times are similar. In the present case we found that this is indeed so. 
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Figure 5: As the interaction well is deepened, i.e. for larger f , unbiased Monte Carlo algorithms 
(top) start to run into trouble. The thermalisation process takes longer and longer as we 
increase the dilution of the system, i.e. by decreasing d. Importance sampling (bottom) 
thermalises much faster: notice the rather large difference in the number of sweeps it takes to 
reach equilibrium for the blue (d = 10 -4 ) curve. 
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MC 
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(N) 






10" 2 


202.4(2) 


7940 


biased 


1CT 3 


199.7(4) 


7302 




10" 4 


193.1(4) 


2857 




lO" 2 


202.0(1) 


41058 


unbiased 


10~ 3 


200.2(1) 


75115 




10~ 4 


196.8(2) 


140683 



Table 1: We estimate the autocorrelation times £ for the instanton number N as a function of 
the Monte Carlo algorithm and the dilution parameter d. The temperature/mass parameter 
is t = 100, i.e. it corresponds to the data in Fig. \5\ We clearly see that as the density of 
the system is lowered, i.e. by decreasing d, the autocorrelation times for unbiased simulations 
become very large compared to the unbiased ones. 

instantons, with roughly equal weight, or settles in an uncorrelated state with 
no pairs. 

For the simulations that follow, we use the 'natural' bonding box, i? max = 
Ro- We always use hot initial conditions, i.e. place the instantons and anti- 
instantons randomly throughout the volume. As we can see from Fig. [4j the 
smallest value of t leads to equivalent equilibrium states for both biased and 
unbiased Monte Carlo. However, as we increase t, and lower d in order to 
maintain a dilute ensemble, ordinary Monte Carlo needs a very long time to 
reach equilibrium compared to biased Monte Carlo, see Fig. [5l Since we start 
with a random distribution, the long thermalisation process can be attributed 
to the fact that it takes longer and longer to locate the attraction centres as 
the system becomes ever more dilute. The correlation between the number of 
instantons N and the number of pairs P is clearly visible in Fig. [6] 

After equilibration, we can estimate autocorrelation times. We find that 
the autocorrelation time in the instanton number iV is much larger for random 
sampling simulations, see Table [TJ The reason is that random insertions and 
deletions are suppressed: insertions because they generally do not form pairs, 
and deletions because they most likely try to break up pairs; this leads to low 
acceptance probabilities. A useful quantity is the autocorrelation time per in- 
stanton. It does not change with the system size, i.e. it is intensive; equivalently, 
autocorrelation times are extensive, i.e. scale linearly with the system size, which 
in our case can be parametrised by the number of instantons. Note that the 
autocorrelation times for the biased simulations actually drop in contrast to 
the unbiased ones. The reason is that as d decreases the system consists of 
a gas of instanton-anti-instanton pairs with less and less interactions, so that 
the updates on pairs will become more and more efficient whereas acceptance 
rates for ordinary Monte Carlo drop ever more@. As the quark mass parameter 



4 Note that in the limit d — ► 0, with t = const, the system behaves as a dilute gas of 
instantons and that ordinary insertions and deletions will work fine again. See Fig. [8]for more 
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Figure 6: The long thcrmalisation time is directly correlated with the slow process of forming 
pairs in the unbiased case. The correlation between pairs and instanton number is also seen 
in the biased simulations; here thcrmalisation is fast (note the different scale, i.e. the upper 
axis) because the algorithm explicitly samples the possible bonding sites. 



t is increased, keeping d fixed, the autocorrelation times increase as well; the 
reason is that the system becomes denser, so that the probability is higher for 
the proposal configurations to be rejected due to repulsive interaction^. 

The toy-model is computationally rather cheap, which allows for very long 
thermalisation sweeps. During such a long equilibration process, helped by the 
low dimensional configuration space, random sampling does eventually form suf- 
ficient pairs to converge to the same state as importance sampling. In higher 
dimensional simulations, as for instance in the IILM, the interaction regions 
have much lower 'entropyH, and it will be harder for ordinary Monte Carlo 
simulations to find the regions of phase space that lead to pair formation. Also, 
realistic systems are computationally more expensiv^], and long runs are pro- 
hibitive. Biased Monte Carlo techniques become unavoidable in these situations. 

The hot initial condition manifests itself by a rapid drop in instanton num- 
bers early in the thermalisation process, due to the lack of pairs; this can even be 
seen in the case of biased simulations, although the algorithm creates pair much 
more quickly, see Fig. [5] We expect that cold initial conditions, i.e. starting 



details. 

5 Note that in the limit t — > oo, with d = const, the system becomes too dense to be treated 
with the methods developed in this paper. 

6 By which we mean that AV/V <C 1, with AV the bonding volume. 

7 For instance, the IfLM in the unquenched case needs to evaluate determinants, a time- 
consuming task. 
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Sweep (biased) 




50 100 150 200 

Sweep (unbiased) 

Figure 7: Random sampling is not well suited to deal with cold initial conditions: we started 
the thermalisation run with a configuration made up of 200 pairs. Random sampling has 
trouble to break up the excess pairs. The inefficiency can be traced back to the fact that 
the algorithm depends solely on the energy difference, which is big when destroying pairs in a 
'naive' way. This leads to low acceptance probabilities and hence long thermalisation sweeps. 
The biased simulations equilibrate very quickly because the algorithm is tailored to deal with 
bonded instantons. 

off with pairs, will also be problematic for random sampling techniques because 
they rely solely on the energy difference between states: in see Fig.[7]we examine 
how ordinary Monte Carlo copes with an ensemble with excess pairs. 

From Figs. [5] and [7] it is clear that biased Monte Carlo techniques vastly 
outperform random sampling. The computational advantage follows from our 
knowledge of the structures that can form, e.g. instanton-anti-instanton molecules 
in the IILM Q , and the ability to construct algorithms that take full advantage 
of that knowledge. 

So far we have restricted the simulations to a single bonding box, which is a 
free parameter a priori. We have run every data set with the two other bonding 
boxes, see Fig. [3] We have found that the results agree on the 2a level. In 
Tabled we show such data for one specific point in parameters space {d, t}. This 
particular set, with large t and very small d 7 is probably representative for more 
realistic simulations in that the unbiased simulations have failed to converge 
to the true equilibrium state; we have lost ergodicity because those regions of 
phase space that correspond to bonded pairs is not sampled adequately. 

The weak dependence on the bonding box allows us to dismiss systematic 
effects, and we can base our choice entirely on achieving low autocorrelation 
times. The data presented in Table [5] would suggest that the two smallest 
bonding boxes are equivalent but the bulk of the computations suggest the 
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Bonding Box (Rq) 


(N) 


(Q 2 ) 


(Sint) 


(0.102328,0.661723) 


194.6(4) 

[1998] 


134(4) 
[1790] 


-2.323(5) 
[1083] 


(0.0457066, 1) 


194.7(4) 
[1907] 


133(3) 
[1212] 


-2.328(7) 
[1812] 


(0.0313807,2) 


195.5(3) 
[4629] 


137(3) 
[2607] 


-2.340(7) 
[5475] 


unbiased 


131.8(3) 
[887] 


129(5) 
[322] 


-0.4(6) 10" b 
[545] 



Table 2: This data is from a simulation with t = 500 and d = 10 -7 , and for these parameters 
the unbiased Monte Carlo simulation equilibrates at a substantially lower value for the total 
number of instantons TV. Note also the low interaction, from which we conclude that no pairs 
have formed. The results for the different bonding boxes agree well at the 2a level. The 
autocorrelations are given in square brackets. The two smallest boxes give rather similar 
autocorrelation times but we have found that in most of the runs the smallest box leads to 
the fastest convergence. 

lowest autocorrelation times are achieved with the smallest bonding box. 

It is interesting to see how the ensemble behaves as we tune the two free 
parameters d and t separately. We have found that the system will always evolve 
to a random state with a negligible amount of pairs upon decreasing d for any 
fixed t. This is as it should be because small d favour a dilute system: from (|28j) 
we see that V increases for decreasing d and constant number of instantons; the 
strong interaction region AV is unaffected by d, and the inequality © will be 
violated for large enough V; thus, random sampling is efficient and the ensemble 
equilibrates in an uncorrelated state. We illustrate this in Fig. [5] 

Similarly, we find that for fixed dilution d and increasingly large interaction 
strength t the system tends to favour energy-dominated configurations. We 
therefore expect that the topological susceptibility should drop to zero since 
pairs do not contribute to the charge fluctuations. The instanton density, on 
the other hand, increases because the energy-dominated configurations favour 
a denser ensemble. This in turn implies that the topological susceptibility will 
be boosted because the topological fluctuations per unit volume have increased. 
Both effects work in opposite directions, but it turns out that the topological 
susceptibility decreases overall, although rather slowly. This is demonstrated in 
Fig. IH1 In particular, the topological susceptibility is lower than the dilute gas 
approximation would suggest. 

In QCD, the level of dilution and the interaction strength change simul- 
taneously, such as to increase both with increasing temperature. Therefore, 
we cannot draw any direct conclusions from this work, other than to note that 
with increasing fermion number the screening, and hence the dilution, decreases 
whereas the fermionic interactions increase. It is therefore not clear a priori 
whether in QCD the IILM will remain in a highly correlated molecular phase 
after the chiral phase transition. 
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Figure 8: The degree of dilution increases from left to right, top to bottom; the temperature- 
mass parameter is fixed at t = 300. As d is decreased fewer instantons are paired up with 
anti-instantons. For the last box the system has equilibrated in a random state. 
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Figure 9: As the interaction strength t is increased, the proportion of pairs to instantons 
increases until basically all instantons are in pairs. Pairs do not contribute to topological 
fluctuations, reducing the topological susceptibility. However, the instanton density increases 
rather strongly with t, resulting in a much denser ensemble. This implies a larger topological 
susceptibility, because the fluctuations per unit volume increase. Both effects are clearly 
antagonistic. In this case, it turns out that the topological susceptibility decreases, although 
rather slowly. 
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5. Conclusions 

We have argued that, quite generically, light quarks in non-trivial back- 
grounds will induce strong and short-ranged interactions at finite temperature. 
The Dirac operator is positive definite and bounded, and leads to an attraction 
between instantons that becomes stronger with decreasing quark mass. Plasma 
effects at finite temperature constrain the spatial extent of the non-trivial back- 
grounds to be below the screening length; since the quark interactions are in- 
duced by the overlap of quark wavcfunctions, centred on top of the classical 
gauge fields, it follows that the interactions become more and more short-ranged 
as the temperature rises. 

Systems with strong and short-ranged interactions are known as strongly 
associating fluids in the literature of computational chemistry and chemical 
engineering. These fluids are hard to simulate and need special techniques. 
Using the concepts of the general purpose Unbonding-Bonding algorithm, we 
extended the biased scheme to grand canonical simulations. Following sugges- 
tions from the literature, we augmented the Monte Carlo steps by pair moves. 
The input from previous studies of the IILM is crucial as these have identified 
the instanton-anti-instanton pairs to be important at finite temperature. This 
knowledge has allowed us to implement biased Monte Carlo techniques that 
specifically address the problems faced with random sampling. 

The defining characteristic of the system is its strong and short-ranged inter- 
actions. We therefore have decided to test importance sampling on a toy-model, 
roughly displaying the features of the IILM. However, we have decided to leave 
out any orientation dependence because it is not necessarily a generic feature 
for the non-trivial backgrounds and would only introduce further complications 
in the Monte Carlo steps. 

The simulations show that random sampling becomes very inefficient if the 
'temperature' is raised or the 'mass' is lowered. We could run the toy-model for 
very long times and check that both the biased and unbiased simulations give 
equivalent results. For one particular set of parameters we actually found that 
ordinary Monte Carlo was not able to reach the correct equilibrium state; in 
this sense we lost ergodicity because the sample was not representative, having 
missed that region of phase space that corresponds to pair formation. We expect 
this to become a much more severe problem for the IILM because the higher 
dimensionality of the latter will make it much harder to sample correctly the 
phase space volume leading to pair formation. Also, the simulations will be 
much more expensive, and long runs will be prohibitive. 

We found that the parameters setting the bonding box have virtually no 
impact on the results, and that bulk properties agree on the 2cr level. This is 
as it should be because a strong dependence on these free parameters would 
introduce some new systematics and render the method impractical. We can 
therefore choose the parameters that give the smallest autocorrelation times. 

Finally, we found that as the degree of dilution is increased, with fixed inter- 
action strength, entropy-dominated configurations will eventually give the bulk 
contribution to the partition function, and the system will settle in a random 
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state. The opposite equilibrium state, namely a highly correlated ensemble of 
pairs, is reached by increasing the interaction strength while holding the dilu- 
tion parameter fixed; in that case the small entropy of these configurations is 
largely compensated by the gain in total energy, and the partition function is 
saturated by these energy-dominated states. These are generic features. In the 
IILM, however, the interactions become stronger and the system more dilute 
simultaneously as the temperature is increased. Whether the system ends up in 
a random ensemble or stays in the molecular phase after the chiral symmetry is 
restored depends on which quantity grows stronger and cannot be decided by 
the results presented in this study. It is, however, interesting to point to the 
possibility of a smaller topological susceptibility in the molecular phase as com- 
pared to the random phase. In light of the ultimate goal of this series of papers, 
this could lead to an axion mass that is rather different from the standard one 
computed in the dilute gas approximation. 
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